SciPy - scientific library in python

Adapted from the notebooks by

Introduction

SciPy builds upon NumPy.

Among others, SciPy deals with:

%matplotlib inline
from scipy import linalg
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from IPython.display import HTML, display

Animation display with python

Animation with matplotlib

You can use FuncAnimation to animate a sequence of images:

%%capture
fig, ax = plt.subplots()
xdata, ydata = [], []
(ln,) = plt.plot([], [], "ro")


def init():
    ax.set_xlim(0, 2 * np.pi)
    ax.set_ylim(-1, 1)
    return (ln,)


def update(frame):
    xdata.append(frame)
    ydata.append(np.sin(frame))
    ln.set_data(xdata, ydata)
    return (ln,)

ani = FuncAnimation(
    fig,
    update,
    interval=50,
    frames=np.linspace(0, 2 * np.pi, 100),
    init_func=init,
    blit=True,
)
display(HTML(ani.to_html5_video()))

Another version exists (using JavaScript):

%%capture
fig, ax = plt.subplots()
xdata, ydata = [], []
(ln,) = plt.plot([], [], "bo")

ani = FuncAnimation(
    fig,
    update,
    interval=50,
    frames=np.linspace(0, 2 * np.pi, 100),
    init_func=init,
    blit=True,
)
display(HTML(ani.to_jshtml()))

More information: http://louistiao.me/posts/notebooks/embedding-matplotlib-animations-in-jupyter-as-interactive-javascript-widgets/

Animation with plotly

matplotlib works fine for advanced tuning, but is harder for simple tasks. So just try plotly for basic animations:

import plotly.express as px
from plotly.offline import plot

df = px.data.gapminder()
fig = px.scatter(
    df,
    x="gdpPercap",
    y="lifeExp",
    animation_frame="year",
    animation_group="country",
    size="pop",
    color="continent",
    hover_name="country",
    log_x=True,
    size_max=55,
    range_x=[100, 100000],
    range_y=[25, 90],
)
fig.show("notebook")

See more here: http://louistiao.me/posts/notebooks/embedding-matplotlib-animations-in-jupyter-as-interactive-javascript-widgets/ ## Linear algebra

scipy for linear algebra : use linalg. It includes functions for solving linear systems, eigenvalues decomposition, SVD, Gaussian elimination (LU, Cholesky), etc.

Documentation: http://docs.scipy.org/doc/scipy/reference/linalg.html

Solving linear systems:

Find x such that: A x = b for specified matrix A and vector b.

 A = np.array([[1, 0, 3], [4, 5, 12], [7, 8, 9]], dtype=float)
# A = np.random.randn(500, 500)
# b = np.ones((500, 1))
b = np.array([[1, 2, 3]], dtype=np.float64).T
# b = np.array([[1, 2, 3]], dtype=np.float64).T

print(A)
print(b)

x = linalg.solve(A, b)
print(x)
print(x.shape)
print(b.shape)
[[ 1.  0.  3.]
 [ 4.  5. 12.]
 [ 7.  8.  9.]]
[[1.]
 [2.]
 [3.]]
[[ 0.8       ]
 [-0.4       ]
 [ 0.06666667]]
(3, 1)
(3, 1)

Check the result at a given precision (different from ==)

np.allclose(A @ x, b, atol=1e-14, rtol=1e-15)
True

Remark: NEVER (or you should really know why) invert a matrix. ALWAYS solve linear systems instead!

Eigen values / vectors

A v_n = \lambda_n v_n with v_n the n-th eigen vector and \lambda_n the n-th eigen value. The associated python functions are eigvals and eig:

A = np.random.randn(3, 3)
A = A + A.T
evals, evecs = linalg.eig(A)
print(evals, "\n ------\n", evecs)

# V = [v_1, ... v_n] , columns = eigen vectors
# A = V diag(s_1,...,s_n) V^T

np.allclose(A, evecs @ np.diag(evals) @ evecs.T)
[ 2.0582851 +0.j -2.66193205+0.j -1.8290659 +0.j] 
 ------
 [[ 0.70513863  0.70663603  0.05869437]
 [ 0.41391987 -0.47742165  0.77507349]
 [-0.57571682  0.5222395   0.62913914]]
True

EXERCISE : Eigen values/vectors

Verify numerically that the output from linalg.eig are indeed approximately eigen values and eigen vectors of matrix A above.

Hint: use https://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html

If A is symmetric you should use eigvalsh (H for Hermitian) instead: this is more robust and leverages the structures (you know they are real!)

Matrix operations

  • linalg.trace(A) # trace
  • linalg.det(A) # determinant
  • linalg.inv(A) # Inverse, consider NEVER using it though :)

Norms

print(linalg.norm(A, ord="fro"))  # fro for Frobenius
print((np.sum(A ** 2)) ** 0.5)
print(linalg.norm(A, ord=2))
print((linalg.eigvalsh(A.T @ A) ** 0.5))
print(linalg.norm(A, ord=np.inf))
3.829869697259075
3.829869697259075
2.661932052341721
[1.8290659  2.0582851  2.66193205]
3.613122169775675

EXERCISE : Norms computation

Check numerically what is the instruction linalg.norm(A, ord=np.inf) computing. Double check with the help and a numerical test.

A = np.random.randn(3, 3)
print(linalg.norm(A, ord=np.inf))
1.9847869788550723

Random generation, distributions, etc.

(see more on this at https://albertcthomas.github.io/good-practices-random-number-generators/ or https://numpy.org/doc/1.18/reference/random/legacy.html#numpy.random.RandomState)

seed = 12345
rng = np.random.default_rng(seed)  # can be called without a seed
rng.random()
0.22733602246716966

Visualization of various popular distributions are available here (see widgets): https://github.com/josephsalmon/Random-Widgets

Optimization

Goal: find functions minima or maxima Documentation: http://scipy-lectures.github.com/advanced/mathematical_optimization/index.html

from scipy import optimize

Finding (local!) minima

def f(x):
    return 4 * x ** 3 + (x - 2) ** 2 + x ** 4


def mf(x):
    return -(4 * x ** 3 + (x - 2) ** 2 + x ** 4)


xs = np.linspace(-5, 3, 100)
plt.figure()
plt.plot(xs, f(xs))
plt.show()

Default solver for minimization/maximization: fmin_bfgs (see https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm)

x_min = optimize.fmin_bfgs(f, x0=-4)
x_max = optimize.fmin_bfgs(mf, x0=-2)
x_min2 = optimize.fmin_bfgs(f, x0=2)


plt.figure()
plt.plot(xs, f(xs))
plt.plot(x_min, f(x_min), "o", markersize=10, color="orange")
plt.plot(x_min2, f(x_min2), "o", markersize=10, color="red")
plt.plot(x_max, f(x_max), "|", markersize=20)
plt.show()
Optimization terminated successfully.
         Current function value: -3.506641
         Iterations: 7
         Function evaluations: 16
         Gradient evaluations: 8
Optimization terminated successfully.
         Current function value: -6.201654
         Iterations: 5
         Function evaluations: 12
         Gradient evaluations: 6
Optimization terminated successfully.
         Current function value: 2.804988
         Iterations: 7
         Function evaluations: 16
         Gradient evaluations: 8

EXERCISE : Bassin of attraction

Draw the points on the curves with two different colors : - orange for the points leading to find the left local minima - red for the points leading to the right local minima.

Find the zeros of a function

Find x such that f(x) = 0, with fsolve.

omega_c = 3.0

def f(omega):
    return np.tan(2 * np.pi * omega) - omega_c / omega



x = np.linspace(1e-8, 3.2, 1000)
y = f(x)

# remove vertical lines when the function flips sign
mask = np.where(np.abs(y) > 50)
x[mask] = y[mask] = np.nan
plt.plot(x, y)
plt.plot([0, 3.3], [0, 0], "k")
plt.ylim(-5, 5)

optimize.fsolve(f, 0.72)
optimize.fsolve(f, 1.1)
optimize.fsolve(f, np.linspace(0.001, 3, 20))
np.unique(np.round(optimize.fsolve(f, np.linspace(0.2, 3, 20)), 3))

my_zeros = (
    np.unique((optimize.fsolve(f, np.linspace(0.2, 3, 20)) * 1000).astype(int)) / 1000.0
)
plt.figure()
plt.plot(x, y, label="$f$")
plt.plot([0, 3.3], [0, 0], "k")
plt.plot(my_zeros, np.zeros(my_zeros.shape), "o", label="$x : f(x)=0$")
plt.legend()
plt.show()

Parameters estimation

from scipy.optimize import curve_fit


def f(x, a, b, c):
    """f(x) = a exp(-bx) + c."""
    return a * np.exp(-b * x) + c


x = np.linspace(0, 4, 50)
y = f(x, 2.5, 1.3, 0.5)  # true signal
yn = y + 0.2 * np.random.randn(len(x))  # noisy added

plt.figure()
plt.plot(x, yn, ".")
plt.plot(x, y, "k", label="$f$")
plt.legend()
plt.show()

(a, b, c), _ = curve_fit(f, x, yn)
print(a, b, c)

2.5880674567536204 1.1582704293784882 0.4083709002776674

Displaying

plt.figure()
plt.plot(x, yn, ".", label="data")
plt.plot(x, y, "k", label="True $f$")
plt.plot(x, f(x, a, b, c), "--k", label="Estimated $\hat{f}$")
plt.legend()
plt.show()

Rem: for polynomial fitting, one can directly use numpy

x = np.linspace(0, 1, 10)
y = np.sin(x * np.pi / 2.0)
line = np.polyfit(x, y, deg=10)
plt.figure()
plt.plot(x, y, ".", label="data")
plt.plot(x, np.polyval(line, x), "k--", label="polynomial approximation")
plt.legend()
plt.show()
/home/jsalmon/anaconda3/envs/peerannot/lib/python3.10/site-packages/IPython/core/interactiveshell.py:3505: RankWarning:

Polyfit may be poorly conditioned

Interpolation

from scipy.interpolate import interp1d, CubicSpline


def f(x):
    return np.sin(x)


n = np.arange(0, 10)
x = np.linspace(0, 9, 100)

y_meas = f(n) + 0.1 * np.random.randn(len(n))  # add noise
y_real = f(x)

linear_interpolation = interp1d(n, y_meas)
y_interp1 = linear_interpolation(x)

cubic_interpolation = CubicSpline(n, y_meas)
y_interp2 = cubic_interpolation(x)


plt.figure()
plt.plot(n, y_meas, "bs", label="noisy data")
plt.plot(x, y_real, "k", lw=2, label="true function")
plt.plot(x, y_interp1, "r", label="linear interp")
plt.plot(x, y_interp2, "g", label="CubicSpline interp")
plt.legend(loc=3)
plt.show()

Images

from scipy import ndimage, misc

img = misc.face()
type(img), img.dtype, img.ndim, img.shape


print(2 ** 8)  # uint8-> code sur 256 niveau.

n_1, n_2, n_3 = img.shape
np.unique(img)


# True image
plt.figure()
plt.imshow(img)
plt.axis("off")
plt.show()
/tmp/ipykernel_25871/4081896913.py:3: DeprecationWarning:

scipy.misc.face has been deprecated in SciPy v1.10.0; and will be completely removed in SciPy v1.12.0. Dataset methods have moved into the scipy.datasets module. Use scipy.datasets.face instead.
256

RGB decomposition

fig, ax = plt.subplots(3, 2)
fig.set_size_inches(9, 6.5)
n_1, n_2, n_3 = img.shape

ax[0, 0].imshow(img[:, :, 0], cmap=plt.cm.Reds)
ax[0, 1].hist(img[:, :, 0].reshape(n_1 * n_2), np.arange(0, 256))

ax[1, 0].imshow(img[:, :, 1], cmap=plt.cm.Greens)
ax[1, 1].hist(img[:, :, 1].reshape(n_1 * n_2), np.arange(0, 256))

ax[2, 0].imshow(img[:, :, 2], cmap=plt.cm.Blues)
ax[2, 1].hist(img[:, :, 2].reshape(n_1 * n_2), np.arange(0, 256))

plt.tight_layout()

print(img.flags)  # cannot edit...
img_compressed = img.copy()
img_compressed.setflags(write=1)
print(img_compressed.flags)  # can edit now


img_compressed[img_compressed < 70] = 50
img_compressed[(img_compressed >= 70) & (img_compressed < 110)] = 100
img_compressed[(img_compressed >= 110) & (img_compressed < 180)] = 150
img_compressed[(img_compressed >= 180)] = 200
plt.figure()
plt.imshow(img_compressed, cmap=plt.cm.gray)
plt.axis("off")
plt.show()
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : False
  WRITEABLE : False
  ALIGNED : True
  WRITEBACKIFCOPY : False

  C_CONTIGUOUS : True
  F_CONTIGUOUS : False
  OWNDATA : True
  WRITEABLE : True
  ALIGNED : True
  WRITEBACKIFCOPY : False

Convert a color image in grayscale

plt.figure()
plt.imshow(np.mean(img, axis=2), cmap=plt.cm.gray)
plt.show()

Blurring (fr: floutage)

img_flou = ndimage.gaussian_filter(img, sigma=20)
fix, ax = plt.subplots()
ax.imshow(img_flou, cmap=plt.cm.gray)
ax.axis("off")
plt.show()

Widget on blurring bandwidth (cannot be rendered online for the moment, only locally)

%matplotlib widget
import ipywidgets as widgets

# set up plot
img_flou = ndimage.gaussian_filter(img, sigma=20)
fix, ax = plt.subplots()
ax.axis("off")
plt.show()

@widgets.interact(sigma=(0.1, 200, 0.1))
def update(sigma=2):
    """Remove old lines from plot and plot new one"""
    # [l.remove() for l in ax.lines]
    img_flou = ndimage.gaussian_filter(img, sigma)
    ax.imshow(img_flou, cmap=plt.cm.gray)

Changing colors in an image

import pooch
import os

url = "https://upload.wikimedia.org/wikipedia/en/thumb/0/05/Flag_of_Brazil.svg/486px-Flag_of_Brazil.svg.png"
name_img =pooch.retrieve(url, known_hash=None)

img = (255 * plt.imread(name_img)).astype(int)
img = img.copy()
plt.figure()
plt.imshow(img[:, :, 2], cmap=plt.cm.gray)


fig, ax = plt.subplots(3, 2)
fig.set_size_inches(9, 6.5)
n_1, n_2, n_3 = img.shape

ax[0, 0].imshow(img[:, :, 0], cmap=plt.cm.Reds)
ax[0, 1].hist(img[:, :, 0].reshape(n_1 * n_2), np.arange(0, 256), density=True)

ax[1, 0].imshow(img[:, :, 1], cmap=plt.cm.Greens)
ax[1, 1].hist(img[:, :, 1].reshape(n_1 * n_2), np.arange(0, 256), density=True)

ax[2, 0].imshow(img[:, :, 2], cmap=plt.cm.Blues)
ax[2, 1].hist(img[:, :, 2].reshape(n_1 * n_2), np.arange(0, 256), density=True)

plt.tight_layout()

EXERCISE : Make the Brazilian italianer

Create a version of the Brazilian flag as follows:

More on colors/ images:

http://josephsalmon.eu/enseignement/Montpellier/HLMA310/matplotlib_slides.pdf

Additional lectures